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ABSTRACT 
In this paper we present an algorithm for adaptively computing 


smooth piecewise polynomial approximations using restricted range 


uniform approximations. We also present several numerical examples and 


offer suggestions for the effective use of this algorithm. We have 
found the algorithm to be effective for approximating a wide class 

of functions, either with or without significant levels of noise. 
Furthermore, since the user of this algorithm actually defines tolerance 


bands within which the approximation will lie, the algorithm allows the 


user a great deal of flexibility and control over the shape of the 


resulting approximations. 
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I. Introduction 
let X be a finite set of real points and let f be a function defined 
on X or let data be given in tabular form. In the 


P P M 
case of data given in tabular form, say {(t., y, te. » we shall set 


1 
X = Coan and define f on X by f(t) 2 es i= iy .ces Me Im what 
follows we shall use this functional notation. Let a = min{x: x € X} 

and b = max{x: x € X}. For any function g defined on X define llell , 

= max{|g(x)|: x € X}. Let SMTH and N be nonnegative integers supplied by 
the user, with N > SMTH. Let u(x) and 2£(x) be functions defined on X 
supplied by the user such that for each x € X we have &(x) < f(x) < u(x) 
and &(x) < u(x). In this setting our algorithm will calculate a piecewise 
polynomial approximation, p, to f, and a set of points iF . © X with 


< «s+ $ % = b such that p restricted to (x; X41) is a 


S| 1 


polynomial p; é Il = {q: q is a real algebraic polynomial of degree < N-1l}, 


N-1 
p has SMTH continuous derivatives (on [a, b]) and for each x EX, 
£(x) < p(x) < u(x). By appropriately choosing £(x) and u(x) then, the 


user obtains an approximation meeting a given preselected (set of) tolerance(s). 


II. The Algorithm 
The algorithm begins by choosing i to be the largest point in X such 
that 
>} 060 Lay x] NX contains at least Nt+l points and 
2) There is a best restricted range uniform approximation Py from 
My) to f (with respect to the constraining curves u(x) and £(x)) 


on [a, x] f. X; that is, there exists PLE My for which 


R(x) < py (x) < u(x) holds for all x © {a, x] fl X and 


lf - p,Il Ca,X,Inx = intl lf - all fa,d, ex’ SS My, 


' £(x) < q(x) < u(x) for all x€ [a,x, 1X}. 


Lf x, = b, then since Py is a piecewise polynomial meeting our 
requirements, we successfully terminate the algorithm. If no such a 
exists, the algorithm fails and an appropriate error message is generated. 
Otherwise, if SMTH = 0 (i.e. we only require the approximation to be 
continuous) we choose the right endpoint of the first subinterval, 


vu 
Xy> to be x If SMTH > 0, in order to add to the stability of the 


1 


by "backing off” from x, in the 


algorithm we (in general) choose x 1 


1 


following manner. 
We first examine the error curve f(x) - P(x) and find those points 


Eye Eye cree by 


choose one of the E's to be X)+ The motivation for choosing x) in this 


in (a, x, IN X at which relative extrema occur. We will 


manner is that in the continuous setting, if f is differentiable and & 

is an interior relative extreme point of f(x) - py, ), then f'(E) - p}(€) = 0 
so that the derivative of Py at &€ would match that of f at €. This 
guarantees that when we smoothly join the next piece of the approximation 

to p, at &, this next piece will closely follow the direction of f at 


ee ; ‘ % 
least near —&. If we merely joined the second piece to the first at x,, no 


2 
such guarantee can be made, and, in fact, severe oscillatory problems tend 

to set in. Our numerical experience indicates that the procedure of 

backing off from x to a smaller Xs contributes very significantly toward 

the stability of the algorithm. To continue, let Pte) be the derivative 

of the centered quadratic interpolation of f evaluated at bye We then choose 


x, to be the largest &, such that [Frce ) - p!(&_)] < EPS, where EPS is a 
v Vv 2 a 


1 
tolerance which can be set by the user, or, if there does not exist such 

a &, then we let x, be the largest §, at which lf"(E,) - py(é,)| is a minimum. 
(In our implementation of the algorithm, we do not in general consider all 


a 
of the relative extreme points of f(x) - p,(x) in [a, x, 1M X, but only 


the largest N - SMTH - 1 of them.) 


We continue by finding successive intervals [x5 xo], [x,, xo], area 

(x,y b] and corresponding polynomial approximations Po» Pgs sees Py €& My _y 
(3) a, ‘ 
to { so that ote 4) x te (x4) for j = 0, 1, ..., SMTH and R(x) <p. (x)<u(x) 
for every x € [x > * Jf\ X for v = 2, ...5) mM, % = b. This is accomplished 
v-1 Vv m 

as follows: 

Suppose we have found the subintervals [a, x], [x), xoJ, sale (x, 5%; 1]; 
and the corresponding approximations Py? Por sees Py_y: Assume further that 
(x, > »] C X contains at least max(2, N-SMTH) points. We now determine an x, and a 


Ps meeting the above requirements. We begin by choosing x; € X to be the 
largest point in X which satisfies 

1) (x; a> x] contains at least max(2, N-SMTH) points and 

2) There exists a best restricted range uniform approximation, Pys 


; z,. fN X, subject to the constraint that Bit tee ) 


1 L 


me ee 
= ppg (Xy_y)9 1 = 05 1, ++» SMTH. 


to £ on fx, 
1- 


x; = b, and the algorithm is successfully 


\ 
If x; = b, we set x; 


: v : ° s ° P 
terminated. If no such x; exists, the algorithm fails and is terminated. 
Otherwise, we choose x; completely analogous to our choice of x) above. 


Finally, we consider the special case where Cx, _ » b] N X contains 


1 


fewer than max(2, N-SMTH) points. Specifically, we choose Sis to bea 


point in X closest to (b - X;_p)/2 which satisfies 


1) Cx 


Xs ys b])M X contains at least max(2, N-SMTH) points and 


2) There exists a best restricted range approximation, Py» to f 


from Myo on (x;y b] 1 X subject to the constraint that 
pe) 5 pyre ~ 5% Oy Ly soos SUTH. 
Again, if we can find such an Xe oy then we change Xia to Xia? set 


x; = b, and successfully terminate the algorithm. If we cannot find such an 


sy the algorithm is terminated and an appropriate error message is generated. 


era i — —_—— 


Remark. In our implementation of this algorithm, the Xs are chosen as 
follows. At each step of this iterative procedure we will let a be 
the current largest point in X such that requirements (1) and (2) 

are satisfied on (x, > 4] X, and we will let Bb be the current 
smallest point in X such that Bb > ¥ and requirement (2) fails to be 
satisfied. We initialize this process by computing (or attempting to 
compute) the best restricted approximation on (x; b] NM X subject to 
the smoothness interpolatory constraints, If this approximation 
satisfies requirement (2), then we set x} = b and we are done. If 

the approximation fails to satisfy (2), we set b= b. Next, let 


t = min{x € X: Cx, _ » x] MX contains at least max(2, N-SMTH) points}. 


1 
If the approximation on [xs _)> t]M X fails to satisfy (2) then the 
algorithm cannot meet the desired accuracy and fails. Otherwise, we 


\ 
set a = t. 


In general, we proceed as follows. We let t = inf{x € X: (b-a)/2<x<b}. 


\ 


If this set is empty, we set t = sup{xe€ X: a <x (S - 3)/2}. it t= a 


fA 


then this procedure has converged and we set x, t =a. Otherwise, 
we compute (or attempt to compute) the best restricted range approxi- 
mation with interpolatory constraints on x;y Pix. Lt this 
approximation satisfies (2) then we set a=t. If this approximation 
fails to satisfy (2) then we set b = t. We continue this process until 
aoe is less than some user definable prescribed tolerance, at which 
point we accept a as a good approximation to A and terminate this 
procedure. We compute the R, in a manner analogous to the above. 

We have implemented a Remes-like algorithm to compute best restricted 


range uniform approximations with interpolatory constraints. See [2] for 


a complete discussion of this problem. 


III. Numerical Results 

By setting u(x) = f(x) + TOL and 2(x) = f(x) - TOL for x € xX, 
where TOL is some positive, preselected tolerance, this algorithm 
simplifies to the best uniform approximation operator version of the 
algorithm given in [3]. For this reason we do not consider here any 
examples with restraining curves differing from the function being 
approximated by a constant. Indeed, the real value of this algorithm 
is that the tolerance we require our approximation to satisfy may vary 
from point to point. Consequently, where f is "nice" we can force our 
approximation to be close to f, and where f is "bad" we can relax our 
requirements, thereby obtaining an approximation which more closely 
reflects the character of f than can be obtained by selecting a fixed 
tolerance throughout the domain of approximation. In the case of 
experimental data which contain considerable levels of noise, the user 
can force the approximation to lie on the "believable" side of the data; 
often, more useful approximations can be obtained with this algorithm 
than can be obtained with (for example) the discrete - version of the 
algorithm given in [3]. 

This algorithm has been implemented as a FORTRAN program running 
on Colorado State University's CDC CYBER 172 and CDC 6400. In the 
appendix we give a listing of the algorithm. As examples, we now 
present approximations to experimental data involving the release of 
bitumen and gas and oil from oil shale heated to a constant temperature 
as a function of time. Because relatively few data points were available 
(14-20), we filled in the gaps between the data points by discretizing 
(200-500 points) the linear interpolation of the original data using an 


algorithm which added (somewhat) more points in regions where the function 


being approximated was complicated (i.e., radically changing slopes 
between data points) and fewer points in regions where the function is 
smooth. The advantage of this unequal spacing over equally spaced points 
is that in regions where the function being approximated is complicated, 
the procedure of “backing off" from a to x; becomes more effective by 
having more densely packed points. Also, the subintervals fe, see) 


may be chosen to be smaller (if needed in order to obtain a close enough 


approximation) while (x, _ f x10 X still contains at least max(2, N-SMTH) 


x 
points. Only the original data points (indicated by "x") are shown in 

the following plots. In each case we chose N = &, and SMTH = 2. The 
curves £(x) and u(x) were chosen by hand or by means of a simple algorithm 
at the original data points, and then they, too, were "filled in" using 
the same linear interpolation scheme as above. The TOL parameter listed 
on the plots is the largest tolerance allowed at any point throughout 


the approximation. This error is not in general reached. 
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Figure 2 


Although the maximum tolerance in both of these examples is fairly large, 
the tolerances throughout most of these intervals of approximation were on the 
order of .15 to 1.5. That is, we required fairly close agreement with the 
function being approximated except at the "bad" points. As an example of what 
can be done by appropriately choosing the restraining curves, we chose a band 
containing the above data and computed restraining curves based on this band 
as the following plot shows. The curve in the center is the resulting 


approximation. 
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Figure 3 


For clarity, we repeated the above plot but without the restraining 


curves. 
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For additional examples using this algorithm and a similar algorithm 
using best L? approximations, see [1]. 

By appropriately setting the restraining curves the user can, 
to a large extent, determine the shape of the approximation. The most 
effective way of determining such restraining curves seems to be trial 
and error. Ideally, these restraining curves would be determined in an 
interactive setting using a graphics terminal with a pen light as follows. 
First one would make a rough initial guess at what the restraining curves 


should be (using some simple algorithm or otherwise), then allow the 


. i eae 


10 


algorithm to compute the first piece of the approximation. One would 
then display the data, the current approximation and the current 
restraining curves and modify the restraining curves on the relevant 
subinterval as desired so that when this first piece is recomputed using 
the updated restraining curves, it behaves as desired. After the user 
is satisfied with the first piece, he would repeat the above strategy 


on each successive subinterval as they are determined by the algorithm. 
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APPENDIX 


Here we give a listing with driver, sample input and output 
which corresponds to Figure 1. The sample input is located after 
the listing of the code and prior to the sample output. This 
example uses the algorithmically defined restraining curves. In 
addition, instructions for user defined restraining curves are 


given in the sample input section. 
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RESTHICTED RANGE ADAFTIVE CURVE FITTING PROGRAM + SAMPLE RUN 
(ALGORITHMICALLY DEFINED RESTRAINING CURVES) 


INPUT 1 
1 ( THIS DENOTES HESTHKAINTS OPTION ) 
6 2 4.00 ol ?5 ( N, SMTH, MAXTOL, MINTOL ) 

( XTABLE, FTABLE) 

3.0 0.0 
5.0 1.3 
720 3.4 TO EMPLOY THE USER DEFINED 
11.0 5.2 RESTRAINING CURVES OPTION, THE 
13.0 6.0 FIRST DATA CARD SHOULD BE 0, 
15.0 14.4 AND THE UPPER AND LOWER TOLER- 
17.5 21.4 ANCES TO BE ALLOWED AT EACH 
20.0 27 4 DATA POINT SHOULD BE ON THE 
22.5 50.9 RESPECTIVE DATA CARDS. 
25.0 49.3 EX: 
eres 47.5 27-5 47.5 140 4.0 
30.0 51.5 (U) (L) 
35.0 36.5 
40.0 27.9 
50.0 94 
60.0 4.2 

OUTPUT 1 


INTEHVAL NUMBER 1 WHICH BEGINS AT .3000000000000000E+01 
AND ENDS AT .1100000000000000F+02 CONTAINS 41 POINTS, 
THE COEFFICIENTS OF BEST APPROXIMATION IN THIS INTERVAL ARE 


C( 1) = .14975652669671 38E+ 02 
C( 2) = -.1144565861352237E+02 
C( 3) = .3049103763268548E+01 
C( 4) = =.3346122135351450E+00 
C( 5) = .1592800176610654E-01 
C( 6) = =.2534785932326346E-03 


THE ERROR OF APPROXIMATION IN THIS INTERVAL IS .3251703402472117E+00. 


INTERVAL NUMBER 2 WHICH BEGINS AT .1100000000000000E+02 
AND ENDS AT .1500000000000000F+02 CONTAINS 24 POINTS. 
THE COEFFICIENTS OF BEST APPROXIMATION IN THIS INTERVAL ARE 


+ 8953993229 386797E+03 
=~ «3853784676856376E+03 
-6639335801 96044 2E+ 02 
- . 5682948501274 041F+01 
+ 2409284901711972E+00 
- «40250421 80238581E-02 


F APPROXIMATION IN THIS INTERVAL IS_~ .1207492591612606E+01. 
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BEST AVAILABLE COPY 
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f 31 
INTERVAL NUMBER 3 WHICH BEGINS AT .1500000000000000E+02 
AND ENDS AT .2000000000000000F+02 CONTAINS 29 POINTS. 
THE COEFFICIENTS OF BEST APPROXIMATION IN THIS INTERVAL ARE 
C »2080819121698872E+05 
= 61199919002 36899h+ O4 
«7147141 380302 528F + 03 
~ 41424944 506629 80E+ 02 
-1192134278058383E+01 


( 
( 
( 
; 
( = .1362678361482650E-01 


) 
) 
) 
) 
) 
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THE ERROR OF APPROXIMATICN IN THIS INTEHKVAL IS .2829293208918557E+01. 


INTERVAL NUMBER 4 WHICH BEGINS AT .2000000000000000E+02 
AND ENDS AT .2359270706355915E+02 CONTAINS 28 POINTS, 
THE COEFFICIENTS OF BEST APPROXIMATION IN THIS INTERVAL ARE 


= «26712586 59287486E+ 06 
-6180344.279850274E+05 
- »5704535104277922E+ 04 
+ 2625284927906778E+ 03 
- .6022618197926448E+01 
55091 32822059026E-01 
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THE ERROR OF APPROXIMATION IN THIS INTERVAL IS-~ .3398666798696013E+01. 


INTERVAL NUMBEK 5 WHICH BEGINS AT .2359270706355915E+02 
AND ENDS AT .2738518411304983E+02 CONTAINS 21 POINTS. 
THE COEFFICIENTS OF BEST APPROXIMATION IN THIS INTERVAL ARE 


-.1013104278693695E+07 
-2002989828887461E+ 06 
- .1582334382597386E+05 
-62437103769 5673 3E+ 03 
-.1230580964377970E+02 
-96914241 8580819 5E-01 
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THE ERROR OF APPROXIMATION IN THIS INTERVAL IS_ .7623136251345386E+00. 


INTERVAL NUMBER 6 WHICH BEGINS AT .2738518411304983E+02 
AND ENDS AT .3343698511741786E+02 CONTAINS 37 POINTS, 
THE COEFFICIENTS OF BEST APPROXIMATION IN THIS INTERVAL ARE 


C( 1) = .2974023810520601E+06 
C( 2) = -.4737982327455143E+05 
C( 2? = .3008544262538955k+04 
C( 4) = ~.9516814694657342E+ 02 
C( 5) = .1499809050599829E+01 
C( 6) = =.9421877586332617E-02 
THE ERHOR OF APPROXIMATION IN THIS INTERVAL IS .1270142945357293E+01. 
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INTEKVAL NUMBER 7 WHICH BEGINS AT .3343698511741786E+02 
AND ENDS AT .40865087038553866E+02 CONTAINS 39 POINTS. 
THE COEFFICIENTS OF BEST APPHOXIMATION IN THIS INTERVAL ARE 


C( 1) = .6567620303936116E+05 
C( 2) = -.8417321059978451E+ 04 
C( 3) = .4305390572534816E+ 03 
C( &) = =.109760057616746 3E+ 02 
C( 5) = .1394356985419112E+00 
C( 6) = -.7061442077258562E-03 


THE ERROR OF APPROXIMATION IN THIS INTERVAL IS” .1270142944891631E+01. 


INTERVAL NUMBER 8 WHICH BEGINS AT .4085087038553866E+02 
AND ENDS AT .6000000000000000E+02 CONTAINS 95 POINTS. 
THE COEFFICIENTS OF BEST APPROXIMATION IN THIS INTERVAL ARE 


C( 1) = .1984768817575497F+05 
C( 2) = -.2082281676169310E+04 
C( 3) = .8723916835701220E+02 
C( 4) = =.1819359280069307E+01 
C( 5) = .1886185858999134E-01 
C( 6) = =.777217759386841 8E- 04 


THE EKRKOR OF APPROXIMATION IN THIS INTERVAL IS .7764477640884024E+00. 


